home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / hermtd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  7.3 KB  |  241 lines
  1. /* linalg/hetd.c
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Factorise a hermitian matrix A into
  21.  *
  22.  * A = U T U'
  23.  *
  24.  * where U is unitary and T is real symmetric tridiagonal.  Only the
  25.  * diagonal and lower triangular part of A is referenced and modified.
  26.  *
  27.  * On exit, T is stored in the diagonal and first subdiagonal of
  28.  * A. Since T is symmetric the upper diagonal is not stored.
  29.  *
  30.  * U is stored as a packed set of Householder transformations in the
  31.  * lower triangular part of the input matrix below the first subdiagonal.
  32.  *
  33.  *  * The full matrix for Q can be obtained as the product
  34.  *
  35.  *       Q = Q_N .. Q_2 Q_1
  36.  *
  37.  * where 
  38.  *
  39.  *       Q_i = (I - tau_i * v_i * v_i')
  40.  *
  41.  * and where v_i is a Householder vector
  42.  *
  43.  *       v_i = [1, A(i+2,i), A(i+3,i), ... , A(N,i)]
  44.  *
  45.  * This storage scheme is the same as in LAPACK.  See LAPACK's
  46.  * chetd2.f for details.
  47.  *
  48.  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
  49.  
  50. #include <config.h>
  51. #include <stdlib.h>
  52. #include <gsl/gsl_math.h>
  53. #include <gsl/gsl_vector.h>
  54. #include <gsl/gsl_matrix.h>
  55. #include <gsl/gsl_blas.h>
  56. #include <gsl/gsl_complex_math.h>
  57.  
  58. #include <gsl/gsl_linalg.h>
  59.  
  60. int 
  61. gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)  
  62. {
  63.   if (A->size1 != A->size2)
  64.     {
  65.       GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
  66.                  GSL_ENOTSQR);
  67.     }
  68.   else if (tau->size + 1 != A->size1)
  69.     {
  70.       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  71.     }
  72.   else
  73.     {
  74.       const size_t N = A->size1;
  75.       size_t i;
  76.   
  77.       const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
  78.       const gsl_complex one = gsl_complex_rect (1.0, 0.0);
  79.       const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);
  80.  
  81.       for (i = 0 ; i < N - 1; i++)
  82.         {
  83.           gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
  84.           gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
  85.           gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);
  86.           
  87.           /* Apply the transformation H^T A H to the remaining columns */
  88.  
  89.           if ((i + 1) < (N - 1) 
  90.               && !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0)) 
  91.             {
  92.               gsl_matrix_complex_view m = 
  93.                 gsl_matrix_complex_submatrix (A, i + 1, i + 1, 
  94.                                               N - (i+1), N - (i+1));
  95.               gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
  96.               gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
  97.               gsl_vector_complex_set (&v.vector, 0, one);
  98.               
  99.               /* x = tau * A * v */
  100.               gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);
  101.  
  102.               /* w = x - (1/2) tau * (x' * v) * v  */
  103.               {
  104.                 gsl_complex xv, txv, alpha;
  105.                 gsl_blas_zdotc(&x.vector, &v.vector, &xv);
  106.                 txv = gsl_complex_mul(tau_i, xv);
  107.                 alpha = gsl_complex_mul_real(txv, -0.5);
  108.                 gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
  109.               }
  110.               
  111.               /* apply the transformation A = A - v w' - w v' */
  112.               gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);
  113.  
  114.               gsl_vector_complex_set (&v.vector, 0, ei);
  115.             }
  116.           
  117.           gsl_vector_complex_set (tau, i, tau_i);
  118.         }
  119.       
  120.       return GSL_SUCCESS;
  121.     }
  122. }  
  123.  
  124.  
  125. /*  Form the orthogonal matrix Q from the packed QR matrix */
  126.  
  127. int
  128. gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A, 
  129.                           const gsl_vector_complex * tau,
  130.                           gsl_matrix_complex * Q, 
  131.                           gsl_vector * diag, 
  132.                           gsl_vector * sdiag)
  133. {
  134.   if (A->size1 !=  A->size2)
  135.     {
  136.       GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
  137.     }
  138.   else if (tau->size + 1 != A->size1)
  139.     {
  140.       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  141.     }
  142.   else if (Q->size1 != A->size1 || Q->size2 != A->size1)
  143.     {
  144.       GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
  145.     }
  146.   else if (diag->size != A->size1)
  147.     {
  148.       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
  149.     }
  150.   else if (sdiag->size + 1 != A->size1)
  151.     {
  152.       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
  153.     }
  154.   else
  155.     {
  156.       const size_t N = A->size1;
  157.  
  158.       size_t i;
  159.  
  160.       /* Initialize Q to the identity */
  161.  
  162.       gsl_matrix_complex_set_identity (Q);
  163.  
  164.       for (i = N - 1; i > 0 && i--;)
  165.     {
  166.           gsl_complex ti = gsl_vector_complex_get (tau, i);
  167.  
  168.           gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (A, i);
  169.  
  170.           gsl_vector_complex_const_view h = 
  171.             gsl_vector_complex_const_subvector (&c.vector, i + 1, N - (i+1));
  172.  
  173.           gsl_matrix_complex_view m = 
  174.             gsl_matrix_complex_submatrix (Q, i + 1, i + 1, N-(i+1), N-(i+1));
  175.  
  176.           gsl_linalg_complex_householder_hm (ti, &h.vector, &m.matrix);
  177.         }
  178.  
  179.       /* Copy diagonal into diag */
  180.  
  181.       for (i = 0; i < N; i++)
  182.         {
  183.           gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
  184.           gsl_vector_set (diag, i, GSL_REAL(Aii));
  185.         }
  186.  
  187.       /* Copy subdiagonal into sdiag */
  188.  
  189.       for (i = 0; i < N - 1; i++)
  190.         {
  191.           gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
  192.           gsl_vector_set (sdiag, i, GSL_REAL(Aji));
  193.         }
  194.  
  195.       return GSL_SUCCESS;
  196.     }
  197. }
  198.  
  199. int
  200. gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A, 
  201.                             gsl_vector * diag, 
  202.                             gsl_vector * sdiag)
  203. {
  204.   if (A->size1 !=  A->size2)
  205.     {
  206.       GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
  207.     }
  208.   else if (diag->size != A->size1)
  209.     {
  210.       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
  211.     }
  212.   else if (sdiag->size + 1 != A->size1)
  213.     {
  214.       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
  215.     }
  216.   else
  217.     {
  218.       const size_t N = A->size1;
  219.  
  220.       size_t i;
  221.  
  222.       /* Copy diagonal into diag */
  223.  
  224.       for (i = 0; i < N; i++)
  225.         {
  226.           gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
  227.           gsl_vector_set (diag, i, GSL_REAL(Aii));
  228.         }
  229.  
  230.       /* Copy subdiagonal into sd */
  231.  
  232.       for (i = 0; i < N - 1; i++)
  233.         {
  234.           gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
  235.           gsl_vector_set (sdiag, i, GSL_REAL(Aji));
  236.         }
  237.  
  238.       return GSL_SUCCESS;
  239.     }
  240. }
  241.